Strong Isotopic Effect in Phase II of Dense Solid Hydrogen and Deuterium 
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Quantum nuclear zero-point motions in solid H2 and D2 under pressure are investigated at 80 K up 
to 160 GPa by first-principles path-integral molecular dynamics calculations. Molecular orientations 
are well-defined in phase II of D2, while solid H2 exhibits large and very asymmetric angular quantum 
fluctuations in this phase, with possible rotation in the (be) plane, making it diflficult to associate 
a well- identified single classical structure. The mechanism for the transition to phase III is also 
described. Existing structural data support this microscopic interpretation. 
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Over the past 20 years, pressure-induced structural 
changes have been studied for most elements up to the 
300 GPa range by a combination of synchrotron diamond 
anvil cell experiments and ab initio calculations. Al- 
though hydrogen has probably been the most extensively 
studied element under pressure, its crystalline structures 
above 10 GPa are still largely unknown. Experimen- 
tal structural determination is hampered by hydrogen's 
weak x-ray scattering cross-section, and the use of neu- 
tron diffraction is equally, if not more, difficult. Theo- 
retically, the difficulties stem from the fact that nuclear 
zero-point motions and exchange contributions can have 
important effects. Novel quantum many-body states of 
matter should be expected, such as the recently predicted 
superconducting/superfluid metallic state [1]. Under- 
standing the phase diagram of hydrogen is thus an open 
challenge in modern condensed matter physics. 

The phase diagrams of solid hydrogen and deuterium 
have been extensively studied by Raman and infra-red 
vibrational spectroscopies. Phase boundaries have been 
determined for three phases The transition from 

phase I to II appears to be related to a quantum order- 
ing of the molecules, whereas that from II to III is more 
likely a classical ordering^. A strong difference in pres- 
sure has been measured between the I-II boundary line 
of solids H2, HD and D2, whereas there is almost none 
for the II-III phase boundary. X-ray diffraction stud- 
ies have shown that the mass centers of the molecules 
remain on an hep lattice at the I-II and II-III transi- 
tions, indicating that these phases are^ differentiated by 
a molecular ordering on this lattice [g*, 7]. The nature of 
this ordering, however, is still unknown. Neutron diffrac- 
tion data may rule out some theoretical predictions for 
phase II of D2 but are insufficient for determination of 
the orientational order [6]. Numerous calculations using 
Density Functional Theory (DFT) have tried to exam- 
ine structural changes in solid molecular hydrogen un- 
der pressure [8]. A very small energy difference between 
various candidate structures is obtained and the relative 
stability of these structures may be easily changed by in- 
cluding zero-point energy contributions. The most recent 
revision of the DFT phase diagram of hydrogen, includ- 
ing zero-point energy at the harmonic level, predicts the 



PGa/m structure for phase II and a C2/c structure for 
phase III [9] though neither of these structures is fully 
compatible with the experimental observations. 

Correct treatment of the nuclear quantum effects in- 
volves calculation of the zero-point motion of the protons 
(deuterons) beyond the harmonic level and to include nu- 
clear exchange to account for the differences between or- 
tho and para molecules. Taking into account the statis- 
tics of the nuclei within an ab initio calculation is still 
a formidable challenge. However, the nuclear quantum 
fluctuations can be accounted for exactly by a calculation 
that combines a DFT treatment of the electrons and a 
path-integral molecular dynamics (PIMD) simulation of 
the nuclear motion. This technique (DFT-PIMD) was 
applied to the case of hydrogen more than a decade ago 
in two studies fiol. It was shown, in both cases, that 
the quasi-harmonic treatment of the zero-point nuclear 
motion is insufficient in high pressure solid hydrogen, and 
the two studies produced different microscopic pictures 
of dense hydrogen. In one case, the quantum fluctuations 
of protons effectively hinder molecular rotation in phase 
II [11] (as a quantum localization) whereas in the other 
case, the structure was identifled as diffuse [loj (coining 
the phrase "quantum fluxional solid" (QFS)). Here we 
revisit a part of the DFT-PIMD phase diagram of hy- 
drogen. Nuclear zero-point motion is shown to have a 
strong influence in phase II: the degree of orientational 
order is important in D2, whereas the structure of solid 
H2 is found consistent with the concept of a QFS. 

We have used the ABINIT code [lH in which we have 
implemented the path-integral formalism for nuclei. In 
this method, the nuclei are treated quantum mechan- 
ically using a discrete representation of the Feynman 
Path-Integral (PI) formulation of quantum statistical me- 
chanics. The PI formalism is based on an isomorphism 
between the quantum system and a classical equivalent 
system. Each quantum nucleus corresponds, in this clas- 
sical equivalent, to a ring polymer of P classical nu- 
clei, with P characterizing the discretization in imagi- 
nary time. The calculations were performed at 80 K. To 
determine the number of imaginary time slices P (a com- 
promise between accuracy and computational burden) , a 
test was performed on H2 at 130 GPa (phase II), using 
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P=64 and P=128. The angular density of probability 
obtained using 128 slices was found to be quasi-identical 
to that obtained using 64, thus P was set to 64 in all the 
calculations. In the strongly harmonic classical equiva- 
lent system, ergodicity is efficiently recovered by using 



tion of state [Hi. A plane -wave cut-off of 20 Hartree was 
required in the present study, in which the pressure was 
limited to 160 GPa and hydrogen remains molecular and 
insulating. 



a Langevin thermostat [IJ 
random number generator 



42 1| including a high-quality 
'22]. For comparison, we have 
also computed classical trajectories by setting P=l. The 
time step is 10 atomic time units (~ 0.24 fs) and 5 in 
the classical case. The mass centers of the molecules 
are initially distributed on an hexagonal lattice with the 
molecules pointing along the c direction. 

It should be noted that our calculation is not a full 
quantum treatment of the problem since we assume dis- 
tinguishable nuclei (Botzmann statistics) and the de- 
coupling between nuclear and electronic degrees of free- 
dom (Born-Oppenheimer approximation). Moreover, we 
emphasize that the objective of the present study is 
not to accurately determine phase transition pressures 
in H2 and D2 (which would require finite-size scaling 
analysis difficult to achieve in ab initio), but rather to 
characterize molecular motions in phase II and its degree 
of orientational order. The DFT-PIMD computations 
were performed in the canonical ensemble, using the ex- 
perimental lattice constants corresponding to pressures 
of 6, 30, 75, 100, 130, 145 and 160 GPa, in a 64-atom su- 
percell for which the Brillouin Zone was sampled by a 2 
X 2 X 2 /c-point mesh (an orthorhombic supercell of size 
2a X 2^/?>a x 2c is constructed, a and c being taken from 
experiments jsl, 0, [l3|)- Imposing the lattice constants 
and a given periodicity is a strong constraint in a system 
where the different molecular structures are very close in 
energy 0. However, even though the structure can be 
impacted, the relative effect - from classical to D2 to H2 
- of quantum fluctuations on the molecular motions, and 
thus the degree of orientational order, should not be de- 
pendent on this constraint, and well reproduced by the 
canonical simulations. 

For the sake of numerical efficiency, we have intro- 
duced an additional level of parallelization in ABINIT 
on the system replicas associated to the discretization 
of the PI in imaginary time, beyond the already exist- 
ing three levels of parallelization (on k-points, bands and 
FFT grids [13]). This level has a quasi-linear scalability 
and, combined with the three others, allows us to per- 
form very long DFT-PIMD trajectories (> 30 000 steps 
with a 64-atom supercell, and in excess of 100 000 steps 
in some cases where a high level of statistics is required) . 
The first 5000 steps are used to equilibrate the system 
and are thus excluded from the averages. The computa- 
tions were performed on the TERA-100 supercomputer 
at the CEA/DAM. The density- functional treatment of 
electrons is achieved in the Generalized Gradient Approx- 
imation (GGA) [14], using the Projector Augmented- 
Wave (PAW) formalism p^. GGA was recently shown 
to give good agreement with experiments for the equa- 
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FIG. 1: Structural properties of D2 under pressure: (1) Pres- 
sure evolution of the ADoP: I (free rotors), II: BSP, III: 
Cmc2i; (2) Pressure evolution of the order parameter; (3) 
ADoP of the classical and D2 structures at 130 GPa, illus- 
trating the effects of quantum fluctuations in D2; (4) View of 
the structure of D2 at 130 GPa (average over 20000 steps) - 
red arrows: primitive cell. 

The DFT-PIMD trajectories are analyzed using dif- 
ferent outputs. We define an order parameter accord- 
ing to [23]: O =< [7^ Eti Er=i ^2(^1z..e,)]' >, with 
P2{X) = ^(3X2 - 1). The brackets < ... > correspond 
to the average over the time steps, N is the number of 
molecules and Qis is the vector of norm 1 lying along 
the i^^ molecule of imaginary time slice (5). The angu- 
lar (normalized) density of probability (ADoP) for the 
(^, (j)) angles of the molecular axis is obtained either by 
averaging over all of the molecules of the supercell, giv- 
ing /(^, 0) or over alternate families of planes along the 
c axis, giving the (normalized) partial angular densities 
of probability /i(<9,(^) (i=l,2). fi{e, (j)) sin {6)606(1) is the 
probality that a molecule of family i points in the (^, (f)) 
direction with uncertainty 60 and ^0. The ADoPs are 
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calculated by averaging over both PIMD time steps and 
imaginary time slices, and therefore include both ther- 
mal and quantum fluctuations. Finally, the structure is 
examined by plotting the atomic positions, averaged over 
the real and imaginary times of the simulation. 

In the classical case, a disordered phase of free rotors 
is found at low pressure while in the pressure range of 
phase II, a complex structure with a 32-atom primitive 
cell is found, identical to the C2/c (idem A2/a) phase 
described by Crespo et al |2J] (see Fig. [l]-4). This phase 
emerges from the Molecular Dynamics, and slightly dif- 
fers from the P6s/m structure found by the most re- 
cent DFT calculations [9] for phase II since it is of lower 
symmetry (Fig. [I]-4). The PGs/m structure has a 16- 
atom unit cell with layers stacked in an ABAB fashion 
whereas the C2/c structure has a 32-atom unit cell with 
layers stacked in the ABA'B' fashion. Both AB and A'B' 
layers show a pinwheel motif. At high pressure 160 
GPa), the C2/c structure of the classical system trans- 
forms into Pca2i. 

The structural evolution of solid D2 is summarized in 
Fig. [H The ADoP at different pressures are compared 
in the first panel. At low pressure, the ADoP is almost 
flat, similar to the one expected from the free rotor state. 
With increasing pressure, clear contours appear in the an- 
gular distribution. Significant changes are seen between 
75 and 100 GPa, corresponding to the I-II transition. In 
addition, strong changes are observed at 145 GPa, cor- 
responding to the II-III transition. The evolution of the 
order parameter with pressure (Fig.[l]-2) also shows these 
two phase transitions. The ADoP of phase II in D2 is 
compared to that of classical hydrogen at 130 GPa in 
Fig. [I1-3. The classical and D2 ADoP look similar, but 
in the case of D2 the ADoP is considerably smoothed by 
the quantum fluctuations. However, despite these strong 
fluctuations, the atomic positions in the supercell yield a 
configuration close to the C2/c structure (Fig.[T]-4). The 
important point is that in this broken-symmetry phase 
(BSP) of D2 at 130 GPa, the molecular orientations are 
well-defined, making it possible to assign a classical struc- 
ture in this pressure range. At 145 GPa, phase II discon- 
tinuously transforms into a Cmc2i structure (phase III) 
with molecules of each basal plane having the same orien- 
tation, alternating along the c axis. All the calculations 
with pressures < 160 GPa exhibit phases in which the 
molecules remain on an hep lattice, in agreement with 
the experiments. 

In the case of hydrogen, the stronger quantum effects 
yield a surprisingly different microscopic picture (Fig.[2|). 
At low pressure, the free rotor phase is associated with 
a flat ADoP. Upon increasing pressure, the molecules 
start to preferentially move in the (be) plane {(j) = 90°), 
as evidenced by the progressive appearance of a broad 
maximum around this value (Fig. [2]-l). Relying only on 
the ADoP /(^,(/)), it is difficult to detect a symmetry 
breaking until 145 GPa. At this pressure, clear changes 




FIG. 2: Structural properties of H2 under pressure: (1) Pres- 
sure evolution of the ADoP: I (free rotors), II: BSP, III: 
Cmc2i; (2) Pressure evolution of the the order parameter; 
(3) QFS at 130 GPa: fi, f2 and integration over of these 
functions, as a function of 0; (4) View of the structure at 130 
GPa (average over 20000 steps). 



are observed: the ADoP exhibits two peaks at Omax ^ 
52.5 and 127.5° (< 6> > = 63 and 117°), with angular 
width (A6>, A0) ^ (28°, 37°): hydrogen evolves towards 
a well-defined Cmc2i structure similar to the one found 
in Ref. 11 at similar pressures. As expected, the quan- 
tum fluctuations in this Cmc2i phase are larger than in 
Cmc2i D2 at the same pressure. On the other hand, 
no clear discontinuity is observed in the evolution of the 
order parameter below this pressure (Fig. [2]-2). 

In fact the symmetry breaking associated with phase II 
occurs at a lower pressure: it is not directly apparent in 
the ADoP /(^, (/)), but in /i and /2 and their (/)-integrated 
counterparts gi{e) and g2{0) (Fig. [13). Below 130 GPa 
these are quasi-identical but at 130 GPa, /i and /2 dif- 
fer significantly, /i has a broad maximum at (^,0) ~ 
(122.5°, 90°) while /2 has a broad maximum at (^, (/>) ~ 
(57.5°, 90°), reflecting a symmetry breaking along the 
c axis that is not evident in the global ADoP /, which 
remains single-peaked. The solid simulated at 130 GPa 
is thus in phase II, while the well-ordered structure at 
145 GPa (Cmc2i) is in phase III. The atomic positions 
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in phase II, obtained by an average over 20000 steps, are 
shown in Fig. [2]-4. No clear structure emerges, reflecting 
the very large quantum fluctuations of protons in phase 
11. For /i and /s, (Al9, A0) ^ (39°, 45°). The fluctu- 
ations in are also very asymmetric, as shown by the 
quantum mechanically averaged values of (^, (/>) ~ (85°, 
90°) and (95°, 90°), which are very different from those 
of the orientation of maximal probability (6 max = 122.5 
and 57.5°). Additional calculations using a 256-atom su- 
percell or P=128 at such pressure provide similar results. 

The averaged structure at 130 GPa is thus close to 
Cmcm - a phase with the molecules lying in the (ab) 
plane - whereas the structure corresponding to the max- 
imum of angular probability would be rather the Cmc2i 
type. This difference is due to large and strongly asym- 
metric angular quantum fluctuations. Two different 
structures may thus be predicted depending on whether 
the most probable molecular orientations or the quantum 
mechanically averaged ones are considered. Moreover, 
inspecting the ADoP of each molecule reveals that they 
can, to some extent, rotate in the (be) plane, showing 
that phase II of H2 cannot be modeled at the harmonic 
level. A few years ago, Biermann et al [10] introduced the 
notion of "quantum fluxional solid" to qualify a solid for 
which an underlying classical structure cannot be clearly 
identified because of large quantum fluctuations. We sug- 
gest this is precisely the case for phase II of hydrogen. 

Upon entering phase III, H2 and D2 adopt the Cmc2i 
structure. At pressures > 160 GPa, the structure in both 
solid isotopes evolve towards a Cmca space group, in 
which the molecular centers are no longer on the hep 
lattice, showing the limit of stability for structures in 
which the molecules remain on the hep sites. However 
a small increase in pressure should result in a displacive 
phase transition. Examination of the nature of phase III 
is precluded at higher pressures as the molecular centers 
move away from an hep configuration and is thus beyond 
the scope of the present work. 
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FIG. 3: Proposed pressure evolution of dense hydrogens at 
T=80 K. 



The proposed phase diagrams at 80 K are summarized 
in Fig. [3] for H2, D2 and the classical case. The most 
important point is that the molecular orientations are 
relatively well-defined in the phase II of D2, which is not 
the case in the phase II of H2, where the concept of QFS 
can be used. There is a strong shift in the I-II transition 
pressure from D2 to H2, in quantitative agreement with 
experiment. In contrast to the DFT-PIMD calculation of 
Ref. [nl, no quantum localization is observed: quantum 
effects systematically induce a gain in symmetry (favor 
molecular rotation). Experimentally, a superstructure 
peak is observed in phase II along the a axis for D2[6j 
but along the c axis for H2[17]. The supercell volume 
was not large enough to permit observation of such super- 
structures in the present simulations. Nevertheless, such 
superstructure peaks indicate that the symmetry break- 
ing is dominant along the c axis for solid H2 and in the 
(ab) plane for solid D2, as suggested here. The very small 
intensity difference in the neutron diffraction peaks mea- 
sured for solid D2 at the I-II transition can be explained 
by the very small difference in the broad ADoP and or- 
der parameter at this transition. Since the change is even 
weaker in the case of solid H2, the structural refinement 
of phase II appears to be beyond the reach of current ex- 
periments. The II-III transition pressure is very similar 
for both isotopes, between 130 and 145 GPa, in qualita- 
tive agreement with experiment (152-165 GPaj^). Fu- 
ture studies will extend the present calculation over the 
pressure domain of phase III in the isothermal-isobaric 
ensemble, thus excluding any assumption regarding the 
lattice constants (unlike in the present work). 
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